library(knitr)
library(tidyverse)
library(haven)
library(RColorBrewer)
library(lubridate)
library(here)
library(zipcode)
library(foreign)
library(lfe)
library(stargazer) 
library(survey)
library(dotwhisker)

fips_kansas_city <- c("Cass" = 29037, "Clay" = 29047, "Jackson" = 29095, "Platte" = 29165)
fips_joplin <- c("Jasper" = 29097, "Newton" = 29145)
fips_nyc <- c("Bronx" = 36005, 
              "Brooklyn/Kings" = 36047, 
              "Manhattan/New York" = 36061, 
              "Queens" = 36081, 
              "Staten Island/Richmond" = 36085)


#####
## Import data and do some basic cleaning of it.
######
## Warshaw-Tausanovitch County Ideology data. I use these to get the state codes for each county. 
county_TW_ideology_estimates<-read.csv(file="county_TW_ideology_estimates.csv")
county_TW_ideology_estimates<-select(county_TW_ideology_estimates, county_fips,county_name, state,abb,mrp_ideology_mean,obama_pc_t2012)
county_TW_ideology_estimates$county_fips<-as.numeric(as.vector(county_TW_ideology_estimates$county_fips))

## Change fips for NYC, Kansas City, and Joplin to the appropriate new value...
county_TW_ideology_estimates$county_fips[county_TW_ideology_estimates$county_fips %in% fips_nyc] <- "36005"
county_TW_ideology_estimates$county_fips[county_TW_ideology_estimates$county_fips %in% fips_kansas_city] <- "Kansas City (MO)"
county_TW_ideology_estimates$county_fips[county_TW_ideology_estimates$county_fips %in% fips_joplin] <- "Joplin (MO)"

states<-group_by(county_TW_ideology_estimates, state, abb)%>%summarise(n=n())
states<-subset(states, state!=".")
states<-select(states,-n)


### Import a file with a list of the dates we'll include in our analysis. We'll use this to expand out the COVID
##  data to generate zeros for days with no COVID deaths/cases
dates<-read.csv(file="dates.csv")
dates<-as.data.frame(dates)
colnames(dates)<-"date"
dates$date<-as.character(dates$date)


### pull number of cases from NYTimes Github
#nyt_county_data_raw <-
#  read_csv(
 #   'https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-counties.csv'
 # ) 
#nyt_county_data_raw$fips<-as.numeric(as.vector(nyt_county_data_raw$fips))
## make backup of NYT data
#write.csv(nyt_county_data_raw, file="data/raw/nyt_county_data_raw.csv")

## For replication purposes, let's use a copy of the data we saved.
nyt_county_data_raw<-read_csv("nyt-covid-us-counties.csv")
nyt_county_data_raw$fips<-as.numeric(as.vector(nyt_county_data_raw$fips))

## Make NYC fips code to the Bronx fip code so it merges properly to survey data 
nyt_county_data_raw$fips[nyt_county_data_raw$county=="New York City"]<-"36005"

# Recode "Kansas City" and Cass, Clay, Jackson, Platte counties to 
# Recode Cass, Clay, Jackson, Platte counties to "Kansas City". Previously,
# Kansas City recoded to Jackson County (29095). This means future data.frames
# will need these changes before merging with NYT data (e.g. Nationscape survey data).
nyt_county_data_raw$fips[nyt_county_data_raw$fips %in% fips_kansas_city] <- "Kansas City (MO)"
nyt_county_data_raw$fips[nyt_county_data_raw$county == "Kansas City"] <- "Kansas City (MO)"
nyt_county_data_raw$fips[nyt_county_data_raw$fips %in% fips_joplin] <- "Joplin (MO)"
nyt_county_data_raw$fips[nyt_county_data_raw$county == "Joplin"] <- "Joplin (MO)"

## Make COVID cases and deaths zeros if they're missing
nyt_county_data_raw$cases[is.na(nyt_county_data_raw$cases)]<-0
nyt_county_data_raw$deaths[is.na(nyt_county_data_raw$deaths)]<-0
nyt_county_data_raw<-select(nyt_county_data_raw, -county)
## Group to ensure each record is distinct at state/fips/date level
nyt_county_data_raw <- nyt_county_data_raw %>% 
  group_by(date, state, fips) %>% 
  summarise(cases = sum(cases),
            deaths = sum(deaths)) %>% 
  ungroup()

# pull county demographics and join to cases by FIPS
#tr_fip_pop_raw <- read_csv('https://raw.githubusercontent.com/tylerreny/county_demos/master/county_pop_1418.csv')
tr_fip_pop_raw  <- read_csv('county_pop_1418.csv')
tr_fip_pop_raw$fips<-as.numeric(as.vector(tr_fip_pop_raw$fips))

## Change fips for NYC, Kansas City, and Joplin to the appropriate new value...
tr_fip_pop_raw$fips[tr_fip_pop_raw$fips %in% fips_nyc] <- "36005"
tr_fip_pop_raw$fips[tr_fip_pop_raw$fips %in% fips_kansas_city] <- "Kansas City (MO)"
tr_fip_pop_raw$fips[tr_fip_pop_raw$fips %in% fips_joplin] <- "Joplin (MO)"
## ... and then aggregate to new fips value level
tr_fip_pop_raw <- tr_fip_pop_raw %>% 
  group_by(fips) %>% 
  summarise(total_pop = sum(total_pop)) %>% 
  ungroup()


state_pop<-merge( tr_fip_pop_raw,county_TW_ideology_estimates, by.x="fips", by.y="county_fips")

state_pop<-group_by(state_pop, state)%>%
	summarise(
		state_total_pop=sum(total_pop, na.rm=T)
		)
dim(state_pop)

## Make a matrix with all dates since beginning of 2020. 
counties_dates<-expand.grid(tr_fip_pop_raw$fips, dates$date)
colnames(counties_dates)<-c("fips", "date")

########
#### Create data for state level Analysis
########

state_analysis<-group_by(nyt_county_data_raw, state,date)%>%
	summarise(	
		state_cases_cumulative=sum(cases, na.rm=T),
		state_deaths_cumulative=sum(deaths, na.rm=T)
		)
dim(state_analysis)

states_dates<-expand.grid(states$state, dates$date)
colnames(states_dates)<-c("state", "date")

state_analysis<-merge(state_analysis, states_dates, by=c("state", "date"), all.y=T)
state_analysis<-merge(state_analysis, states, by="state")
state_analysis$state_cases_cumulative[is.na(state_analysis$state_cases_cumulative)]<-0
state_analysis$state_deaths_cumulative[is.na(state_analysis$state_deaths_cumulative)]<-0
state_analysis<-merge(state_analysis, state_pop, by="state")

state_analysis <- group_by(state_analysis,state) %>%
  mutate(
    state_deaths_10days=state_deaths_cumulative-lag(state_deaths_cumulative,10, default = 0), ## deaths over past 10 days
    state_cases_10days=state_cases_cumulative-lag(state_cases_cumulative,10, default = 0),## cases over past 10 days
    state_deaths_20days=state_deaths_cumulative-lag(state_deaths_cumulative,20, default = 0), ## deaths over past 20 days
    state_cases_20days=state_cases_cumulative-lag(state_cases_cumulative,20, default = 0),## cases over past 20 days
    state_deaths_30days=state_deaths_cumulative-lag(state_deaths_cumulative,30, default = 0), ## deaths over past 30 days
    state_cases_30days=state_cases_cumulative-lag(state_cases_cumulative,30, default = 0),## cases over past 30 days
    state_deaths_40days=state_deaths_cumulative-lag(state_deaths_cumulative,40, default = 0), ## deaths over past 40 days
    state_cases_40days=state_cases_cumulative-lag(state_cases_cumulative,40, default = 0),## cases over past 40 days
    state_deaths_50days=state_deaths_cumulative-lag(state_deaths_cumulative,50, default = 0), ## deaths over past 50 days
    state_cases_50days=state_cases_cumulative-lag(state_cases_cumulative,50, default = 0),## cases over past 50 days
    state_deaths_60days=state_deaths_cumulative-lag(state_deaths_cumulative,60, default = 0), ## deaths over past 60 days
    state_cases_60days=state_cases_cumulative-lag(state_cases_cumulative,60, default = 0),## cases over past 60 days
    state_deaths_70days=state_deaths_cumulative-lag(state_deaths_cumulative,70, default = 0), ## deaths over past 70 days
    state_cases_70days=state_cases_cumulative-lag(state_cases_cumulative,70, default = 0),## cases over past 70 days
    state_deaths_80days=state_deaths_cumulative-lag(state_deaths_cumulative,80, default = 0), ## deaths over past 80 days
    state_cases_80days=state_cases_cumulative-lag(state_cases_cumulative,80, default = 0),## cases over past 80 days
    state_deaths_90days=state_deaths_cumulative-lag(state_deaths_cumulative,90, default = 0), ## deaths over past 90 days
    state_cases_90days=state_cases_cumulative-lag(state_cases_cumulative,90, default = 0),## cases over past 90 days
    state_deaths_30days_lead=lead(state_deaths_cumulative,30)-state_deaths_cumulative, ## deaths over next 30 days
    state_deaths_90days_lead=lead(state_deaths_cumulative,90)-state_deaths_cumulative ## deaths over next 90 days
  )

state_analysis$state_deaths_10days[state_analysis$state_deaths_10days<0]<-0
state_analysis$state_deaths_20days[state_analysis$state_deaths_20days<0]<-0
state_analysis$state_deaths_30days[state_analysis$state_deaths_30days<0]<-0
state_analysis$state_deaths_40days[state_analysis$state_deaths_40days<0]<-0
state_analysis$state_deaths_50days[state_analysis$state_deaths_50days<0]<-0
state_analysis$state_deaths_60days[state_analysis$state_deaths_60days<0]<-0
state_analysis$state_deaths_70days[state_analysis$state_deaths_70days<0]<-0
state_analysis$state_deaths_80days[state_analysis$state_deaths_80days<0]<-0
state_analysis$state_deaths_90days[state_analysis$state_deaths_90days<0]<-0
state_analysis$state_deaths_30days_lead[state_analysis$state_deaths_30days_lead<0]<-0
state_analysis$state_deaths_90days_lead[state_analysis$state_deaths_90days_lead<0]<-0

state_analysis <- group_by(state_analysis,state) %>%
	mutate(	
		state_deaths_per_100000_10days=state_deaths_10days/state_total_pop*100000,
		state_deaths_per_100000_20days=state_deaths_20days/state_total_pop*100000,
		state_deaths_per_100000_30days=state_deaths_30days/state_total_pop*100000,
		state_deaths_per_100000_40days=state_deaths_40days/state_total_pop*100000,
		state_deaths_per_100000_50days=state_deaths_50days/state_total_pop*100000,
		state_deaths_per_100000_60days=state_deaths_60days/state_total_pop*100000,
		state_deaths_per_100000_70days=state_deaths_70days/state_total_pop*100000,
		state_deaths_per_100000_80days=state_deaths_80days/state_total_pop*100000,
		state_deaths_per_100000_90days=state_deaths_90days/state_total_pop*100000,
		state_deaths_per_100000_30days_lead=state_deaths_90days/state_total_pop*100000,
		state_deaths_per_100000_90days_lead=state_deaths_90days_lead/state_total_pop*100000,
		state_deaths_per_100000_log_10days=log(state_deaths_per_100000_10days+.1),
		state_deaths_per_100000_log_20days=log(state_deaths_per_100000_20days+.1),
		state_deaths_per_100000_log_30days=log(state_deaths_per_100000_30days+.1),
		state_deaths_per_100000_log_40days=log(state_deaths_per_100000_40days+.1),
		state_deaths_per_100000_log_50days=log(state_deaths_per_100000_50days+.1),
		state_deaths_per_100000_log_60days=log(state_deaths_per_100000_60days+.1),
		state_deaths_per_100000_log_70days=log(state_deaths_per_100000_70days+.1),
		state_deaths_per_100000_log_80days=log(state_deaths_per_100000_80days+.1),
		state_deaths_per_100000_log_90days=log(state_deaths_per_100000_90days+.1),
		state_deaths_per_100000_log_30days_lead=log(state_deaths_per_100000_30days_lead+.1),
		state_deaths_per_100000_log_90days_lead=log(state_deaths_per_100000_90days_lead+.1),
		state_deaths_per_100000_log_cumulative=log(state_deaths_cumulative/state_total_pop*100000+.1)
		
		)


nyt_county_data_raw<-select(nyt_county_data_raw, -state)
nyt_county_data_raw<-merge(nyt_county_data_raw, counties_dates, by=c("fips", "date"), all.y=T)

# Impute missing cases and deaths as 0s. That is, if the fips/date combination
# does not appear in the NYT dataset, code it as no cases/deaths occurring.
nyt_county_data_raw$cases[is.na(nyt_county_data_raw$cases)]<-0
nyt_county_data_raw$deaths[is.na(nyt_county_data_raw$deaths)]<-0

########
#### Create data for state level Analysis
########

county_analysis<-merge( nyt_county_data_raw,county_TW_ideology_estimates, by.x="fips", by.y="county_fips")

county_analysis<-merge(county_analysis, tr_fip_pop_raw, by="fips")
county_analysis<-group_by(county_analysis, fips,date)%>%
	summarise(	
		county_total_pop=sum(total_pop, na.rm=T),
		county_cases_cumulative=sum(cases, na.rm=T),
		county_deaths_cumulative=sum(deaths, na.rm=T)
		)
dim(county_analysis)

county_analysis <- group_by(county_analysis,fips) %>%
	mutate(
		county_deaths_10days=county_deaths_cumulative-lag(county_deaths_cumulative,10, default = 0), ## deaths over past 10 days
		county_cases_10days=county_cases_cumulative-lag(county_cases_cumulative,10, default = 0),## cases over past 10 days
		county_deaths_20days=county_deaths_cumulative-lag(county_deaths_cumulative,20, default = 0), ## deaths over past 20 days
		county_cases_20days=county_cases_cumulative-lag(county_cases_cumulative,20, default = 0),## cases over past 20 days
		county_deaths_30days=county_deaths_cumulative-lag(county_deaths_cumulative,30, default = 0), ## deaths over past 30 days
		county_cases_30days=county_cases_cumulative-lag(county_cases_cumulative,30, default = 0),## cases over past 30 days
		county_deaths_40days=county_deaths_cumulative-lag(county_deaths_cumulative,40, default = 0), ## deaths over past 40 days
		county_cases_40days=county_cases_cumulative-lag(county_cases_cumulative,40, default = 0),## cases over past 40 days
		county_deaths_50days=county_deaths_cumulative-lag(county_deaths_cumulative,50, default = 0), ## deaths over past 50 days
		county_cases_50days=county_cases_cumulative-lag(county_cases_cumulative,50, default = 0),## cases over past 50 days
		county_deaths_60days=county_deaths_cumulative-lag(county_deaths_cumulative,60, default = 0), ## deaths over past 60 days
		county_cases_60days=county_cases_cumulative-lag(county_cases_cumulative,60, default = 0),## cases over past 60 days
		county_deaths_70days=county_deaths_cumulative-lag(county_deaths_cumulative,70, default = 0), ## deaths over past 70 days
		county_cases_70days=county_cases_cumulative-lag(county_cases_cumulative,70, default = 0),## cases over past 70 days
		county_deaths_80days=county_deaths_cumulative-lag(county_deaths_cumulative,80, default = 0), ## deaths over past 80 days
		county_cases_80days=county_cases_cumulative-lag(county_cases_cumulative,80, default = 0),## cases over past 80 days
		county_deaths_90days=county_deaths_cumulative-lag(county_deaths_cumulative,90, default = 0), ## deaths over past 90 days
		county_cases_90days=county_cases_cumulative-lag(county_cases_cumulative,90, default = 0)## cases over past 90 days
	)


county_analysis$county_deaths_10days[county_analysis$county_deaths_10days<0]<-0
county_analysis$county_deaths_20days[county_analysis$county_deaths_20days<0]<-0
county_analysis$county_deaths_30days[county_analysis$county_deaths_30days<0]<-0
county_analysis$county_deaths_40days[county_analysis$county_deaths_40days<0]<-0
county_analysis$county_deaths_50days[county_analysis$county_deaths_50days<0]<-0
county_analysis$county_deaths_60days[county_analysis$county_deaths_60days<0]<-0
county_analysis$county_deaths_70days[county_analysis$county_deaths_70days<0]<-0
county_analysis$county_deaths_80days[county_analysis$county_deaths_80days<0]<-0
county_analysis$county_deaths_90days[county_analysis$county_deaths_90days<0]<-0
county_analysis <- group_by(county_analysis,fips) %>%
	mutate(	
		county_deaths_per_100000_10days=county_deaths_10days/county_total_pop*100000,
		county_deaths_per_100000_20days=county_deaths_20days/county_total_pop*100000,
		county_deaths_per_100000_30days=county_deaths_30days/county_total_pop*100000,
		county_deaths_per_100000_40days=county_deaths_40days/county_total_pop*100000,
		county_deaths_per_100000_50days=county_deaths_50days/county_total_pop*100000,
		county_deaths_per_100000_60days=county_deaths_60days/county_total_pop*100000,
		county_deaths_per_100000_70days=county_deaths_70days/county_total_pop*100000,
		county_deaths_per_100000_80days=county_deaths_80days/county_total_pop*100000,
		county_deaths_per_100000_90days=county_deaths_90days/county_total_pop*100000,
		county_deaths_per_100000_log_10days=log(county_deaths_per_100000_10days+.1),
		county_deaths_per_100000_log_20days=log(county_deaths_per_100000_20days+.1),
		county_deaths_per_100000_log_30days=log(county_deaths_per_100000_30days+.1),
		county_deaths_per_100000_log_40days=log(county_deaths_per_100000_40days+.1),
		county_deaths_per_100000_log_50days=log(county_deaths_per_100000_50days+.1),
		county_deaths_per_100000_log_60days=log(county_deaths_per_100000_60days+.1),
		county_deaths_per_100000_log_70days=log(county_deaths_per_100000_70days+.1),
		county_deaths_per_100000_log_80days=log(county_deaths_per_100000_80days+.1),
		county_deaths_per_100000_log_90days=log(county_deaths_per_100000_90days+.1)
		)

## Import a cleaned, trimmed version of the Nationscape survey data
survey<-read.dta( "nationscape_covid_analysis.dta")

# Impute fips to match NYT + population data
survey$county_backup <- survey$county
survey$county<-as.numeric(as.vector(survey$county))
survey$county[survey$county %in% fips_nyc] <- "36005"
survey$county[survey$county %in% fips_kansas_city] <- "Kansas City (MO)"
survey$county[survey$county %in% fips_joplin] <- "Joplin (MO)"

## create a variable for Biden vs. Trump voting intention
survey$vote_pres<-NA
survey$vote_pres[survey$trump_biden==1]<-0
survey$vote_pres[survey$trump_biden==2]<-1

survey$vote_pres_dk<-NA
survey$vote_pres_dk[survey$trump_biden==1]<-0 # Biden
survey$vote_pres_dk[survey$trump_biden==2]<-1 # Trump
survey$vote_pres_dk[survey$trump_biden==999]<-.5

survey$vote_house<-NA
survey$vote_house[survey$house_intent==1]<-0 # Democrat
survey$vote_house[survey$house_intent==2]<-1 # Republican

survey$vote_house_dk<-NA
survey$vote_house_dk[survey$house_intent==1]<-0
survey$vote_house_dk[survey$house_intent==2]<-1
survey$vote_house_dk[survey$house_intent==999]<-.5

## make senate vote variables and make states without senate races NA
survey$vote_senate<-NA
survey$vote_senate[survey$senate_intent==1]<-0 # Democrat
survey$vote_senate[survey$senate_intent==2]<-1 # Republican
survey$vote_senate[
survey$state=="CA" | 
survey$state=="CT" |
survey$state=="DC" |
survey$state=="FL" |
survey$state=="HI" |
survey$state=="IN" | 
survey$state=="MD" |
survey$state=="MO" | 
survey$state=="ND" | 
survey$state=="NV" | 
survey$state=="NY" |
survey$state=="OH" |
survey$state=="PA" |
survey$state=="UT" | 
survey$state=="VT" |
survey$state=="WA" | 
survey$state=="WI"  ]<-NA

survey$vote_senate_dk<-NA
survey$vote_senate_dk[survey$senate_intent==1]<-0
survey$vote_senate_dk[survey$senate_intent==2]<-1
survey$vote_senate_dk[survey$senate_intent==999]<-.5
survey$vote_senate_dk[
survey$state=="CA" | 
survey$state=="CT" |
survey$state=="DC" |
survey$state=="FL" |
survey$state=="HI" |
survey$state=="IN" | 
survey$state=="MD" |
survey$state=="MO" | 
survey$state=="ND" | 
survey$state=="NV" | 
survey$state=="NY" |
survey$state=="OH" |
survey$state=="PA" |
survey$state=="UT" | 
survey$state=="VT" |
survey$state=="WA" | 
survey$state=="WI" ]<-NA

# Took out regions<-read.csv("bea_regions.csv"). Not used.

### Create simplified presidential approval variables
survey$pres_approval2<-NA
survey$pres_approval2[survey$pres_approval==1| survey$pres_approval==2]<-1
survey$pres_approval2[survey$pres_approval==3| survey$pres_approval==4]<-0

survey$pres_approval2_dk<-NA
survey$pres_approval2_dk[survey$pres_approval==1| survey$pres_approval==2]<-1
survey$pres_approval2_dk[survey$pres_approval==3| survey$pres_approval==4]<-0
survey$pres_approval2_dk[survey$pres_approval==999]<-.5


## Merge the data on COVID at the state level with our survey data
dim(survey)
survey$date<-as.character(survey$date)
state_analysis$date<-as.character(state_analysis$date)
dim(survey)
state_analysis$abb<-as.vector(state_analysis$abb)
survey<-merge(survey, state_analysis, by.x=c("state", "date"),by.y=c("abb", "date"), all.x=T)
dim(survey)

survey$county<-as.character(survey$county)
county_analysis$fips<-as.character(county_analysis$fips)
county_analysis$date<-as.character(county_analysis$date)
survey<-merge(survey,county_analysis, by.x=c("county", "date"),by.y=c("fips", "date"), all.x=T)

## Make presidential approval on a 0-100 scale to make it easier to interpret
survey$pres_approval2<-survey$pres_approval2*100
survey$pres_approval2_dk<-survey$pres_approval2_dk*100
survey$vote_pres<-survey$vote_pres*100
survey$vote_pres_dk<-survey$vote_pres_dk*100
survey$vote_house<-survey$vote_house*100
survey$vote_house_dk<-survey$vote_house_dk*100
survey$vote_senate<-survey$vote_senate*100
survey$vote_senate_dk<-survey$vote_senate_dk*100

survey$cluster_county<-paste(survey$date, survey$county, sep="-")
survey$cluster_state<-paste(survey$date, survey$state, sep="-")


## descriptive stats
do_viz<-0
if(do_viz==1){
viz<-group_by(survey, state, month, year)%>%
summarise(
	state_total_pop=mean(state_total_pop, na.rm=T),
	state_deaths_cumulative=mean(state_deaths_cumulative, na.rm=T),
	state_deaths_per_100000_log_30days=mean(state_deaths_per_100000_log_30days, na.rm=T),
	state_deaths_per_100000_log_60days=mean(state_deaths_per_100000_log_60days, na.rm=T),
	state_deaths_per_100000_log_cumulative=mean(state_deaths_per_100000_log_cumulative, na.rm=T)	
	)
viz<-arrange(viz, year, month)
as.data.frame(viz[viz$county==36005,])
}

# -------------------------------------------------------------------- #

survey$race4<-NA
survey$race4[survey$race_ethnicity==1]<-"white"
survey$race4[survey$race_ethnicity==2]<-"black"
survey$race4[survey$race_ethnicity>2]<-"other"
survey$race4[survey$hispanic!=1]<-"hispanic"

survey$education5<-NA
survey$education5[survey$education<4]<-"no hs"
survey$education5[survey$education==4]<-"hs"
survey$education5[survey$education>4]<-"some college"
survey$education5[survey$education>7]<-"college"
survey$education5[survey$education>8]<-"grad school"

survey$race_education_gender_vote2016<-paste(survey$race4, survey$education5, survey$gender, survey$vote_2016, sep="-")

## Create index variables
survey<-mutate(survey,
	republican_preference=rowMeans(cbind(pres_approval2,vote_pres, vote_senate, vote_house), na.rm=T),
	republican_preference_dk=rowMeans(cbind(pres_approval2_dk,vote_pres_dk, vote_senate_dk, vote_house_dk), na.rm=T)
	)

#####
## Regressions
######

summary(survey$state_deaths_per_100000_30days)
summary(survey$county_deaths_per_100000_30days)
summary(survey$state_deaths_per_100000_log_30days)
summary(survey$county_deaths_per_100000_log_30days)

######
### Republican preference
######

### State and county, binary outcome variable
reg_index_state_log<-(felm(republican_preference~(state_deaths_per_100000_log_30days)|wave+state+race_education_gender_vote2016|0|cluster_state, data=survey, weight=survey$weight))
reg_index_county_log<-(felm(republican_preference~(county_deaths_per_100000_log_30days)|wave+county+race_education_gender_vote2016|0|cluster_county, data=survey, weight=survey$weight))

######
### presidential approval as outcome
######

### state Level, binary outcome variable
reg_approval2_state_log<-(felm(pres_approval2~(state_deaths_per_100000_log_30days)|wave+state+race_education_gender_vote2016|0|cluster_state, data=survey, weight=survey$weight))

### county Level, binary outcome variable
reg_approval2_county_log<-(felm(pres_approval2~(county_deaths_per_100000_log_30days)|wave+county+race_education_gender_vote2016|0|cluster_county, data=survey, weight=survey$weight))

######
### Presidential vote as outcome
######

### state Level, binary outcome variable
reg_votepres_state_log<-(felm(vote_pres~(state_deaths_per_100000_log_30days)|wave+state+race_education_gender_vote2016|0|cluster_state, data=survey, weight=survey$weight))

### county Level, binary outcome variable
reg_votepres_county_log<-(felm(vote_pres~(county_deaths_per_100000_log_30days)|wave+county+race_education_gender_vote2016|0|cluster_county, data=survey, weight=survey$weight))

######
### House vote as outcome
######

### state Level, binary outcome variable
reg_votehouse_state_log<-(felm(vote_house~(state_deaths_per_100000_log_30days)|wave+state+race_education_gender_vote2016|0|cluster_state, data=survey, weight=survey$weight))

### county Level, binary outcome variable
reg_votehouse_county_log<-(felm(vote_house~(county_deaths_per_100000_log_30days)|wave+county+race_education_gender_vote2016|0|cluster_county, data=survey, weight=survey$weight))

######
### Senate vote as outcome
######

### state Level, binary outcome variable
reg_votesenate_state_log<-(felm(vote_senate~(state_deaths_per_100000_log_30days)|wave+state+race_education_gender_vote2016|0|cluster_state, data=survey, weight=survey$weight))

### county Level, binary outcome variable
reg_votesenate_county_log<-(felm(vote_senate~(county_deaths_per_100000_log_30days)|wave+county+race_education_gender_vote2016|0|cluster_county, data=survey, weight=survey$weight))

####
## Out a dot and whiskers plot of our main results
####

## with index
#stargazer(reg_approval2_county_log,reg_approval2_state_log,reg_votepres_county_log,reg_votepres_state_log,reg_votesenate_county_log,reg_votesenate_state_log,reg_votehouse_county_log,reg_votehouse_state_log,reg_index_county_log,reg_index_state_log)
#results<-matrix(data=NA, nrow=10, ncol=0)
#results<-as.data.frame(results)
#results$term<-c("Trump Approval", "Trump Approval", "President Vote", "President Vote", "Senate Vote", "Senate Vote", "House Vote", "House Vote", "Index", "Index")
#results$estimate<-c(coef(reg_approval2_county_log)[1],coef(reg_approval2_state_log)[1],coef(reg_votepres_county_log)[1],coef(reg_votepres_state_log)[1],coef(reg_votesenate_county_log)[1],coef(reg_votesenate_state_log)[1],coef(reg_votehouse_county_log)[1],coef(reg_votehouse_state_log)[1],coef(reg_index_county_log)[1],coef(reg_index_state_log)[1])
#results$std.error<-c(summary(reg_approval2_county_log)$coefficients[2],summary(reg_approval2_state_log)$coefficients[2],summary(reg_votepres_county_log)$coefficients[2],summary(reg_votepres_state_log)$coefficients[2],summary(reg_votesenate_county_log)$coefficients[2],summary(reg_votesenate_state_log)$coefficients[2],summary(reg_votehouse_county_log)$coefficients[2],summary(reg_votehouse_state_log)$coefficients[2],summary(reg_index_county_log)$coefficients[2],summary(reg_index_state_log)$coefficients[2])
#results$model<-c("County", "State", "County", "State", "County", "State","County", "State","County", "State")

## without index
stargazer(reg_approval2_county_log,reg_approval2_state_log,reg_votepres_county_log,reg_votepres_state_log,reg_votesenate_county_log,reg_votesenate_state_log,reg_votehouse_county_log,reg_votehouse_state_log)
results<-matrix(data=NA, nrow=8, ncol=0)
results<-as.data.frame(results)
results$term<-c("Trump Approval", "Trump Approval", "President Vote", "President Vote", "Senate Vote", "Senate Vote", "House Vote", "House Vote")
results$estimate<-c(coef(reg_approval2_county_log)[1],coef(reg_approval2_state_log)[1],coef(reg_votepres_county_log)[1],coef(reg_votepres_state_log)[1],coef(reg_votesenate_county_log)[1],coef(reg_votesenate_state_log)[1],coef(reg_votehouse_county_log)[1],coef(reg_votehouse_state_log)[1])
results$std.error<-c(summary(reg_approval2_county_log)$coefficients[2],summary(reg_approval2_state_log)$coefficients[2],summary(reg_votepres_county_log)$coefficients[2],summary(reg_votepres_state_log)$coefficients[2],summary(reg_votesenate_county_log)$coefficients[2],summary(reg_votesenate_state_log)$coefficients[2],summary(reg_votehouse_county_log)$coefficients[2],summary(reg_votehouse_state_log)$coefficients[2])
results$model<-c("County", "State", "County", "State", "County", "State","County", "State")


## Convert units to a doubling of fatalities instead of one unit on log-scale
results$estimate<-results$estimate*.693
results$std.error<-results$std.error*.693
results$upper<-round(results$estimate+results$std.error*1.96,2)
results$lower<-round(results$estimate-results$std.error*1.96,2)


pdf("Figures/figure2_main_results_short.pdf", height=4)
dwplot(results)+ theme_bw()+ theme(legend.position="bottom")  + theme(legend.title = element_blank())+ xlab("Effect of doubling in COVID deaths per 100,000 people") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=12))+ theme(legend.text=element_text(size=14))
dev.off()

#####
### Robustness
#####


####
## Testing different clustering strategies
####

## states
robustness_clustering<-matrix(data=NA, nrow=3, ncol=4)
robustness_clustering<-as.data.frame(robustness_clustering)
colnames(robustness_clustering)<-c("term", "estimate", "std.error", "p.value")

temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+state+race_education_gender_vote2016|0|state, data=subset(survey), weight=survey$weight))
robustness_clustering$term[1]<-"States"
robustness_clustering$estimate[1]<-coef(temp)[1]
robustness_clustering$std.error[1]<-summary(temp)$coefficients[2]
robustness_clustering$p.value[1]<-summary(temp)$coefficients[4]
	
temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
robustness_clustering$term[2]<-"State-Day"
robustness_clustering$estimate[2]<-coef(temp)[1]
robustness_clustering$std.error[2]<-summary(temp)$coefficients[2]
robustness_clustering$p.value[2]<-summary(temp)$coefficients[4]

temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+state+race_education_gender_vote2016|0|0, data=subset(survey), weight=survey$weight))
robustness_clustering$term[3]<-"None"
robustness_clustering$estimate[3]<-coef(temp)[1]
robustness_clustering$std.error[3]<-summary(temp)$coefficients[2]
robustness_clustering$p.value[3]<-summary(temp)$coefficients[4]

pdf("Figures/robustness_clustering_states.pdf", height=2.5)
dwplot(robustness_clustering)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=8)) + theme(legend.position="none")
dev.off()

## counties
robustness_clustering<-matrix(data=NA, nrow=3, ncol=4)
robustness_clustering<-as.data.frame(robustness_clustering)
colnames(robustness_clustering)<-c("term", "estimate", "std.error", "p.value")

temp<-(felm(republican_preference~county_deaths_per_100000_log_30days|wave+county+race_education_gender_vote2016|0|county, data=subset(survey), weight=survey$weight))
robustness_clustering$term[1]<-"States"
robustness_clustering$estimate[1]<-coef(temp)[1]
robustness_clustering$std.error[1]<-summary(temp)$coefficients[2]
robustness_clustering$p.value[1]<-summary(temp)$coefficients[4]

temp<-(felm(republican_preference~county_deaths_per_100000_log_30days|wave+county+race_education_gender_vote2016|0|cluster_county, data=subset(survey), weight=survey$weight))
robustness_clustering$term[2]<-"State-Day"
robustness_clustering$estimate[2]<-coef(temp)[1]
robustness_clustering$std.error[2]<-summary(temp)$coefficients[2]
robustness_clustering$p.value[2]<-summary(temp)$coefficients[4]

temp<-(felm(republican_preference~county_deaths_per_100000_log_30days|wave+county+race_education_gender_vote2016|0|0, data=subset(survey), weight=survey$weight))
robustness_clustering$term[3]<-"None"
robustness_clustering$estimate[3]<-coef(temp)[1]
robustness_clustering$std.error[3]<-summary(temp)$coefficients[2]
robustness_clustering$p.value[3]<-summary(temp)$coefficients[4]

pdf("Figures/robustness_clustering_counties.pdf", height=3)
dwplot(robustness_clustering)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=5)) + theme(legend.position="none")
dev.off()

####
## No control variables
####
robustness_FE<-matrix(data=NA, nrow=3, ncol=4)
robustness_FE<-as.data.frame(robustness_FE)
colnames(robustness_FE)<-c("term", "estimate", "std.error", "p.value")

temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
robustness_FE$term[1]<-"Main Specification: 
Time + State + Individual"
robustness_FE$estimate[1]<-coef(temp)[1]
robustness_FE$std.error[1]<-summary(temp)$coefficients[2]
robustness_FE$p.value[1]<-summary(temp)$coefficients[4]

temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+state|0|cluster_state, data=subset(survey), weight=survey$weight))
robustness_FE$term[2]<-"Time + State"
robustness_FE$estimate[2]<-coef(temp)[1]
robustness_FE$std.error[2]<-summary(temp)$coefficients[2]
robustness_FE$p.value[2]<-summary(temp)$coefficients[4]

temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
robustness_FE$term[3]<-"Time + Individual"
robustness_FE$estimate[3]<-coef(temp)[1]
robustness_FE$std.error[3]<-summary(temp)$coefficients[2]
robustness_FE$p.value[3]<-summary(temp)$coefficients[4]

pdf("Figures/robustness_FEs_state.pdf", height=2.3)
dwplot(robustness_FE)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=8)) + theme(legend.position="none")
dev.off()

main_county<-(felm(republican_preference~county_deaths_per_100000_log_30days|wave+county+race_education_gender_vote2016|0|cluster_county, data=subset(survey), weight=survey$weight))
nocontrols_county<-(felm(republican_preference~county_deaths_per_100000_log_30days|wave+county|0|cluster_county, data=subset(survey), weight=survey$weight))
nocountyFE_county<-(felm(republican_preference~county_deaths_per_100000_log_30days|wave+race_education_gender_vote2016|0|cluster_county, data=subset(survey), weight=survey$weight))
stargazer(main_county, nocontrols_county,nocountyFE_county)


####
## Dropping states
####
robustness_states<-matrix(data=NA, nrow=50, ncol=4)
robustness_states<-as.data.frame(robustness_states)
colnames(robustness_states)<-c("term", "estimate", "std.error", "p.value")

for(i in 1:50){
	print(i)
	temp<-(felm(republican_preference~(state_deaths_per_100000_log_30days)|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey, state!=states$abb[i]), weight=survey$weight[survey$state!=states$abb[i]]))
	robustness_states$term[i]<-states$state[i]
	robustness_states$estimate[i]<-coef(temp)[1]
	robustness_states$std.error[i]<-summary(temp)$coefficients[2]
	robustness_states$p.value[i]<-summary(temp)$coefficients[4]
}
robustness_states$term<-states$state
pdf("Figures/robustness_dropping_states_state.pdf", height=5.5)
dwplot(robustness_states)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=5)) + theme(legend.position="none")
dev.off()

###
## Different lags for treatment
####
## State-level
treatments<-c("state_deaths_per_100000_log_10days", 
	"state_deaths_per_100000_log_20days",
	"state_deaths_per_100000_log_30days",
	"state_deaths_per_100000_log_40days",
	"state_deaths_per_100000_log_50days",
	"state_deaths_per_100000_log_60days",
	"state_deaths_per_100000_log_70days",
	"state_deaths_per_100000_log_80days",
	"state_deaths_per_100000_log_90days")
treatments2<-c("10 Days", "20 Days", "30 Days", "40 Days",
"50 Days", "60 Days", "70 Days", "80 Days","90 Days", "Cumulative")

robustness_lags<-matrix(data=NA, nrow=length(treatments2), ncol=4)
robustness_lags<-matrix(data=NA, nrow=length(treatments2), ncol=4)
robustness_lags<-as.data.frame(robustness_lags)
colnames(robustness_lags)<-c("term", "estimate", "std.error", "p.value")
	
for(i in 1:length(treatments)){
	print(i)
	treatment<-treatments[i]
	temp<-(felm(republican_preference~get(treatment)|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
	robustness_lags$term[i]<-treatments[i]
	robustness_lags$estimate[i]<-coef(temp)[1]
	robustness_lags$std.error[i]<-summary(temp)$coefficients[2]
	robustness_lags$p.value[i]<-summary(temp)$coefficients[4]
}
robustness_lags$term<-treatments2
## cumulative
temp<-(felm(republican_preference~state_deaths_per_100000_log_cumulative|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
	robustness_lags$term[10]<-treatments2[10]
	robustness_lags$estimate[10]<-coef(temp)[1]
	robustness_lags$std.error[10]<-summary(temp)$coefficients[2]
	robustness_lags$p.value[10]<-summary(temp)$coefficients[4]

pdf("Figures/robustness_lags_state.pdf", height=2.75)
dwplot(robustness_lags)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=8)) + theme(legend.position="none")
dev.off()

## County-level
treatments<-c("county_deaths_per_100000_log_10days", 
	"county_deaths_per_100000_log_20days",
	"county_deaths_per_100000_log_30days",
	"county_deaths_per_100000_log_40days",
	"county_deaths_per_100000_log_50days",
	"county_deaths_per_100000_log_60days",
	"county_deaths_per_100000_log_70days",
	"county_deaths_per_100000_log_80days",
	"county_deaths_per_100000_log_90days")
treatments2<-c("10 Days", "20 Days", "30 Days", "40 Days",
"50 Days", "60 Days", "70 Days", "80 Days","90 Days")
robustness_lags<-matrix(data=NA, nrow=length(treatments), ncol=4)
robustness_lags<-as.data.frame(robustness_lags)
colnames(robustness_lags)<-c("term", "estimate", "std.error", "p.value")
	
for(i in 1:length(treatments)){
	print(i)
	treatment<-treatments[i]
	temp<-(felm(republican_preference~get(treatment)|wave+county+race_education_gender_vote2016|0|cluster_county, data=subset(survey), weight=survey$weight))
	robustness_lags$term[i]<-treatments2[i]
	robustness_lags$estimate[i]<-coef(temp)[1]
	robustness_lags$std.error[i]<-summary(temp)$coefficients[2]
	robustness_lags$p.value[i]<-summary(temp)$coefficients[4]
}

robustness_lags$term<-treatments2

pdf("Figures/robustness_lags_county.pdf", height=4.5)
dwplot(robustness_lags)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=8)) + theme(legend.position="none")
dev.off()


###
## Different outcomes
####
## State-level
outcomes<-c("republican_preference", 
	"republican_preference_dk","pres_approval2", "pres_approval2_dk", "vote_pres", "vote_pres_dk", "vote_senate", "vote_senate_dk", "vote_house", "vote_house_dk" )
outcomes2<-c("Index", 
	"Index","Pres. Approval", "Pres. Approval", "Pres. Vote", "Pres. Vote", "Senate Vote", "Senate Vote", "House Vote", "House Vote" )
models2<-c("No DKs", "With DKs","No DKs", "With DKs","No DKs", "With DKs","No DKs", "With DKs","No DKs", "With DKs")

robustness_state_dks<-matrix(data=NA, nrow=length(outcomes2), ncol=5)
robustness_state_dks<-matrix(data=NA, nrow=length(outcomes2), ncol=5)
robustness_state_dks<-as.data.frame(robustness_state_dks)
colnames(robustness_state_dks)<-c("model","term", "estimate", "std.error", "p.value")
	
for(i in 1:length(models2)){
	print(i)
	outcome<-outcomes[i]
	temp<-(felm(get(outcome)~state_deaths_per_100000_log_30days|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
	robustness_state_dks$model[i]<-models2[i]
	robustness_state_dks$term[i]<-outcomes2[i]
	robustness_state_dks$estimate[i]<-coef(temp)[1]
	robustness_state_dks$std.error[i]<-summary(temp)$coefficients[2]
	robustness_state_dks$p.value[i]<-summary(temp)$coefficients[4]
}

pdf("Figures/robustness_state_dks.pdf", height=2.75)
dwplot(robustness_state_dks)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=8)) 
dev.off()


## County-level
outcomes<-c("republican_preference", 
	"republican_preference_dk","pres_approval2", "pres_approval2_dk", "vote_pres", "vote_pres_dk", "vote_senate", "vote_senate_dk", "vote_house", "vote_house_dk" )
outcomes2<-c("Index", 
	"Index","Pres. Approval", "Pres. Approval", "Pres. Vote", "Pres. Vote", "Senate Vote", "Senate Vote", "House Vote", "House Vote" )
models2<-c("No DKs", "With DKs","No DKs", "With DKs","No DKs", "With DKs","No DKs", "With DKs","No DKs", "With DKs")

robustness_county_dks<-matrix(data=NA, nrow=length(outcomes2), ncol=5)
robustness_county_dks<-matrix(data=NA, nrow=length(outcomes2), ncol=5)
robustness_county_dks<-as.data.frame(robustness_county_dks)
colnames(robustness_county_dks)<-c("model","term", "estimate", "std.error", "p.value")
	
for(i in 1:length(models2)){
	print(i)
	outcome<-outcomes[i]
	temp<-(felm(get(outcome)~county_deaths_per_100000_log_30days|wave+county+race_education_gender_vote2016|0|cluster_county, data=subset(survey), weight=survey$weight))
	robustness_county_dks$model[i]<-models2[i]
	robustness_county_dks$term[i]<-outcomes2[i]
	robustness_county_dks$estimate[i]<-coef(temp)[1]
	robustness_county_dks$std.error[i]<-summary(temp)$coefficients[2]
	robustness_county_dks$p.value[i]<-summary(temp)$coefficients[4]
}

pdf("Figures/robustness_county_dks.pdf", height=2.75)
dwplot(robustness_county_dks)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=8)) 
dev.off()



### Placebo checks

placebo<-matrix(data=NA, nrow=3, ncol=4)
placebo<-matrix(data=NA, nrow=3, ncol=4)
placebo<-as.data.frame(placebo)
colnames(placebo)<-c("term", "estimate", "std.error", "p.value")
	
temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
	placebo$term[1]<-"Main Specification
	(deaths over past 30 days)"
	placebo$estimate[1]<-coef(temp)[1]
	placebo$std.error[1]<-summary(temp)$coefficients[2]
	placebo$p.value[1]<-summary(temp)$coefficients[4]

temp<-(felm(republican_preference~ state_deaths_per_100000_log_30days_lead |wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey, ( month<4) | year==2019), weight=survey$weight[( survey$month<4) | survey $year==2019]))
	placebo$term[2]<-"30 Day Lead"
	placebo$estimate[2]<-coef(temp)[1]
	placebo$std.error[2]<-summary(temp)$coefficients[2]
	placebo$p.value[2]<-summary(temp)$coefficients[4]
temp<-(felm(republican_preference~ state_deaths_per_100000_log_90days_lead |wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey, ( month<4) | year==2019), weight=survey$weight[( survey$month<4) | survey $year==2019]))
	placebo$term[3]<-"90 Day Lead"
	placebo$estimate[3]<-coef(temp)[1]
	placebo$std.error[3]<-summary(temp)$coefficients[2]
	placebo$p.value[3]<-summary(temp)$coefficients[4]


pdf("Figures/placebo.pdf", height=2.75)
dwplot(placebo)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0,.25,.5,.75,1), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%",".25%", ".5%",".75%", "1%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=12),axis.text=element_text(size=8)) + theme(legend.position="none")
dev.off()



###
## Cross sectional model
###

robustness_XS<-matrix(data=NA, nrow=2, ncol=4)
robustness_XS<-as.data.frame(robustness_XS)
colnames(robustness_XS)<-c("term", "estimate", "std.error", "p.value")

temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))

temp<-(felm(republican_preference~state_deaths_per_100000_log_cumulative|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
robustness_XS$term[1]<-"Main Specification using cumulative deaths"
robustness_XS$estimate[1]<-coef(temp)[1]
robustness_XS$std.error[1]<-summary(temp)$coefficients[2]
robustness_XS$p.value[1]<-summary(temp)$coefficients[4]

temp<-(felm(republican_preference~state_deaths_per_100000_log_cumulative|wave+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
robustness_XS$term[2]<-"No state FE using cumulative deaths"
robustness_XS$estimate[2]<-coef(temp)[1]
robustness_XS$std.error[2]<-summary(temp)$coefficients[2]
robustness_XS$p.value[2]<-summary(temp)$coefficients[4]

pdf("Figures/robustness_xs.pdf", height=2.5)
dwplot(robustness_XS)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=10),axis.text=element_text(size=8)) + theme(legend.position="none")
dev.off()



###
## Comparison of cumulative deaths and main results
###

robustness_war<-matrix(data=NA, nrow=2, ncol=4)
robustness_war<-as.data.frame(robustness_war)
colnames(robustness_war)<-c("term", "estimate", "std.error", "p.value")

temp<-(felm(republican_preference~state_deaths_per_100000_log_30days|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
robustness_war$term[1]<-"Main Specification, with 30 day cutoff for deaths"
robustness_war$estimate[1]<-coef(temp)[1]
robustness_war$std.error[1]<-summary(temp)$coefficients[2]
robustness_war$p.value[1]<-summary(temp)$coefficients[4]

temp<-(felm(republican_preference~state_deaths_per_100000_log_cumulative|wave+state+race_education_gender_vote2016|0|cluster_state, data=subset(survey), weight=survey$weight))
robustness_war$term[2]<-"Main Specification using cumulative deaths"
robustness_war$estimate[2]<-coef(temp)[1]
robustness_war$std.error[2]<-summary(temp)$coefficients[2]
robustness_war$p.value[2]<-summary(temp)$coefficients[4]


pdf("Figures/robustness_war.pdf", height=2.5)
dwplot(robustness_war)+ theme_bw() + xlab("Effect of 1-unit increase in log(COVID deaths per 100,000 people)") +scale_x_continuous( breaks=c(-1.75,-1.5,-1.25,-1,-.75,-.5,-.25,0), labels=c("-1.75%","-1.5%","-1.25%","-1%","-.75%","-.5%","-.25%","0%"))+geom_vline(xintercept=0,linetype="dotted")+theme(axis.title=element_text(size=10),axis.text=element_text(size=8)) + theme(legend.position="none")
dev.off()
